%% 载波频率

% 为简化计算，我们产生了PRBS7的信号，其占12700个点即127位每位100个点，其采样频率为10e6 Hz。
% 而具体的实验中，我们使用的是1064nm的种子，其频率约为3*10e14 Hz。而实际信号的频率为800Mbps和1Gbps，即8*10e8和10e9
% 则实际情况下，一个PRBS信号频率为800/127=6.29M Hz，时长为1.59*e-7 s。
% 所以如果使用一个完整的PRBS7的周期作为计算，则光的频率对应为3*10e14 Hz  /6.29M =4.77e7  角频率为2 \pi *  4.77e7
% =3e8。
%% 信号频率

fs=10^6;%采样频率
T=1/fs;%采样点数
L=12700;%信号长度
t=(0:L-1)*T;%时间序列 
f=fs*(0:(L-1))/L;
binwidth = fs/length(t);
f_hz=[0:binwidth:fs-binwidth];%将点数转化为频率
%% 定义载波

w=3e8;%载波角频率
Z1=exp(1i*w*t);%载波
plot(t(1:500),Z1(1:500))
plot(f_hz,abs(fftshift(fft(Z1))));
%% 
%% 读取PRBS信号

load("Z.mat");%读取之前处理好的PRBS信号
plot(t,Z)
% plot(f_hz,abs(fftshift(fft(Z))));
plot(f_hz,abs(fft(Z)));
title('Amplitude Spectrum of PRBS7')
xlabel('f (Hz)')
ylabel('|P1(f)|')
PRBS_F=fft(Z);
plot(f_hz,PRBS_F)

a=0.13;%用来调制PRBS信号的高频分量。
PRBS_F(a*length(PRBS_F):(1-a)*length(PRBS_F))=0; %将频率中的高频部分变为0
plot(f_hz,PRBS_F);
%% 画出傅里叶谱的一半
P2 = PRBS_F/L;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(L/2))/L;
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
PRBS_lb=ifft(PRBS_F);
plot(t,PRBS_lb);
title('PRBS7')
xlabel('t')
ylabel('Signal strength')
%% 调制

% y=exp(1i.*(2.*pi.*f_c.*t+beta.*pm))
beta=pi;%调制深度
Z2=exp(1i*w*t+1i*beta*Z)+randn(size(t))/10;
plot(t,Z2);
fft_z2=fft(Z2);
plot(f_hz,abs(fftshift(fft_z2))/(length(Z2)/2));
title('PRBS7-FFT')
xlabel('f(Hz)')
ylabel('Signal strength')
Z_RC=exp(1i*w*t+1i*beta*PRBS_lb)+randn(size(t))/10;
% plot(t,Z_RC);
fft_z2_RC=fft(Z_RC);
 plot(f_hz,abs(fftshift(fft_z2_RC))/(length(fft_z2_RC)/2));
 title('PRBS7-FFT')
xlabel('f(Hz)')
ylabel('Signal strength')
 FFTC=abs(fftshift(fft_z2_RC))/(length(fft_z2_RC)/2);
 %% 将图像平移到中心，方便看图
C=max(FFTC,[],[1,2]);%获得极大值点
B=find(FFTC==C);%寻找极大值点对应的横坐标
FFTC1=FFTC(1:B*2-1);
tll=(0:B*2-1-1)*T;
f_hz_tll=[0:binwidth:tll-binwidth]
plot(tll,FFTC1);
plot(tll,20*log10(FFTC1));
title('PRBS7-FFT')
xlabel('f')
ylabel('Signal strength(dB)')
%% 倍频

Z22=exp(2*(1i*w*t+1i*beta*Z))+randn(size(t))/10;
plot(t,Z22);
fft_z22=fft(Z22);
plot(t,abs(fftshift(fft_z22))/(length(Z22)));
Z12=exp(2*(1i*w*t))+randn(size(t))/10;
plot(t,Z12);
fft_z12=fft(Z12);
plot(t,abs(fftshift(fft_z12))/(length(Z12)/2));

Z_RC2=exp(2*(1i*w*t+1i*beta*PRBS_lb))+randn(size(t))/10;
plot(t,Z_RC2);
fft_z2_RC2=fft(Z_RC2);
plot(f_hz,abs(fftshift(fft_z2_RC2))/(length(fft_z2_RC2)/2));
title('SHG-PRBS7-FFT')
xlabel('f(Hz)')
ylabel('Signal strength')
%  FFTC2=abs(fftshift(fft_z2_RC2))/(length(fft_z2_RC2)/2);
% C2=max(FFTC2,[],[1,2]);%获得极大值点
% B2=find(FFTC2==C2);%寻找极大值点对应的横坐标
% FFTC12=FFTC2(1:B2*2-1);
% tll2=(0:B2*2-1-1)*T;
% plot(tll2,FFTC12);
% plot(tll2,20*log10(FFTC12));
%% 添加sin
% sin函数

N1=64;
si=0.5.*sin(2*pi*N1*t/(T*L))+0.5;
si2=0.5.*sin(2*pi*N1*t/(T*L))+0.5;
plot(t,si);
% beta=pi/2;
% Z1_sin=exp(1i*w*t+2*beta*1i*si);
% plot(t,Z1_sin);
% plot(f_hz,abs(fftshift(fft(Z1_sin))));
% Z_RC_sin=exp(1i*w*t+1i*beta*PRBS_lb+1i*2*beta*si);%+randn(size(t))/10;
Z_RC_sin=exp(1i*w*t+1i*beta*si);
plot(t,Z_RC_sin);
fft_z2_sin=fft(Z_RC_sin);
plot(f_hz,abs(fftshift(fft_z2_sin))/(length(Z_RC_sin)/2));
title('Sin-FFT')
xlabel('f(Hz)')
ylabel('Signal strength')

%% 方波

N1=64;
pm=0.5.*square(2*pi*N1*t/(T*L))+0.5;
plot(t,pm);
% Z_RC_pm=exp(1i*w*t+1i*beta*PRBS_lb+1i*beta*pm);
Z_RC_pm=exp(1i*w*t+1i*beta*pm)+randn(size(t))/10;
plot(t,Z_RC_pm);
fft_z2_pm=fft(Z_RC_pm);
plot(f_hz,abs(fftshift(fft_z2_pm))/(length(Z_RC_pm)/2));
title('Square-FFT')
xlabel('f(Hz)')
ylabel('Signal strength')
% plot(pm);%产生了N1个方波 用来对比PRBS信号
%% 三角波

san= 0.5*sawtooth(2*pi*N1*t/(T*L))+0.5;
plot(t,san);
san_fft=fft(san);
plot(t,abs(san_fft));
Z_san=exp(1i*w*t+1i*beta*san);
plot(t,Z_san);
Z_san_fft=fft(Z_san);
plot(t,abs(Z_san_fft)/(length(Z_san)/2));
title('Sawtooth-FFT')
xlabel('f(Hz)')
ylabel('Signal strength')